
import sys
import indexing.featureindex
import coverage.basecountreader

if __name__=="__main__":
    sample1 = sys.argv[1]
    sample2 = sys.argv[2]
    if len(sys.argv)<=3:
        path1 = "/archive10/bowtie/Sample_Parameters/%s/Current" % sample1
        path2 = "/archive10/bowtie/Sample_Parameters/%s/Current" % sample2
    else:
        path1 = sys.argv[3]
        path2 = sys.argv[4]
    coverageindex1 = "%s/%s.coverage_summary.INDEX" % (path1, sample1)
    coverageindex2 = "%s/%s.coverage_summary.INDEX" % (path2, sample2)
    ix1 = indexing.featureindex.FeatureIndex(coverageindex1)
    ix2 = indexing.featureindex.FeatureIndex(coverageindex2)
    reader1 = coverage.basecountreader.BaseCountReader(ix1)
    reader2 = coverage.basecountreader.BaseCountReader(ix2)

    n=0
    cts = {}
    for f1 in ix1.getFeatures():
        (chr, pos) = (f1.chr, f1.spos)
        f2 = reader2[(chr, pos)]
        if f2:
            f1 = f1.getDetails()
            f2 = f2.getDetails()
            key = (min(f1.reads, 100), min(f2.reads, 100))
            cts[key] = cts.get(key, 0) + 1
        n+=1
        if 0==(n%100000):
            print >> sys.stderr, "# ", n, "of", ix1.getFeatureCount(), ":", len(cts)
    keys = cts.keys()
    keys.sort()
    for (n1, n2) in keys:
        print "%d\t%d\t%ld" % (n1, n2, cts[(n1, n2)])
